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Abstract 

Flows through a transonic diffuser were investigated 
with the PARC code using five turbulence models to 
determine the effects of turbulence model selection on flow 
prediction. Three of the turbulence models were algebraic 
models: Thomas (the standard algebraic turbulence model 
in PARC), Baldwin-Lomax, and Modified Mixing Length- 
Thomas (MMLT). The other two models were the low 
Reynolds number k-e models of Chien and Speziale. Three 
diffuser flows, referred to as the no-shock, weak-shock, and 
strong-shock cases, were calculated with each model to 
conduct the evaluation. Pressure distributions, velocity 
profiles, locations of shocks, and maximum Mach numbers 
in the duct were the flow quantities compared. Overall, the 
Chien k-e model was the most accurate of the five models 
when considering results obtained for all three cases. 
However, the MMLT model provided solutions as accurate 
as the Chien model for the no-shock and the weak-shock 
cases, at a substantially lower computational cost (measured 
in CPU time required to obtain converged solutions). The 
strong shock flow, which included a region of shock- 
induced flow separation, was only predicted well by the two 
k-8 models. 

Nomenclature 

H diffuser height (varying through duct) 

HLj. diffuser height at throat 

k turbulent kinetic energy 

P static pressure 

P o inflow total pressure (also reference pressure) 

P o standard deviation of pressure 

R ratio of outflow static pressure to inflow total 
pressure, P/P o 

U ref reference velocity (speed of sound based on 
reference temperature and pressure) 

V velocity 

x,y Cartesian coordinates 

y* nondimensional vertical position 

£ rate of turbulent kinetic energy dissipation 


Introduction 

Computational fluid dynamics (CFD) is now being 
used extensively to analyze flows through advanced 
propulsion systems. These flows often include characteris- 
tics such as attached and separated turbulent boundary 
layers, oblique and normal shocks, shock wave/boundary 
layer interactions, turbulent mixing, and other complex 
phenomena. The most sophisticated CFD codes employing 
Navier-Stokes solvers are required to analyze propulsion 
components with flow characteristics such as these. Despite 
the advances in flow solving capabilities, the ability of 
Navier-Stokes solvers to calculate complex flows is strongly 
dependent on the turbulence model employed. 

In the current study, flows through a transonic diffuser 
were calculated with the PARC code, a general purpose 
Navier-Stokes solver for fluid flow simulation, using five of 
the turbulence models installed in PARC. Calculations 
obtained with PARC using each turbulence model were 
compared with experimental data to determine the effects of 
turbulence model selection on the prediction of diffuser 
flows. The flow quantities under comparison were the 
pressure distributions along the top and bottom walls of the 
diffuser, velocity profiles, locations of shocks in the flows, 
and Mach numbers in the duct. The computational cost 
required to obtain solutions using the different turbulence 
models was also considered. 

The following sections describe the PARC code and 
turbulence models used in the study, diffuser cases that were 
examined, and comparison of PARC calculations with 
experimental data. 


The PARC Code 

The PARC code 1 - 2 is an internal flow Navier-Stokes 
code used extensively by government and industry to 
analyze propulsion flows, especially those of aircraft engine 
inlets and nozzles. PARC was derived from the ARC 
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external flow Navier-Stokes code. 3 - 4 One version of the 
PARC code contains the two-dimensional and axisymmetric 
solver (PARC2D) whereas the other version contains the 
three-dimensional solver (PARC3D). The governing 
equations of motion are the time-dependent Reynolds 
averaged Navier-Stokes equations satisfying a perfect gas 
relationship and Fourier's heat conduction law. These 
equations are discretized in conservation law form with 
respect to general curvilinear coordinates and solved with 
the Beam and Warming approximate factorization algo- 
rithm. 5 Although a time-dependent solver based on the 
work of Jameson 6 is available, PARC is intended for steady- 
state flow simulations. 

Turbulence Models in PARC 

Three of the turbulence models investigated in this 
study are algebraic (zero-equation) models and two are 
energy-dissipation rate (k-£) two-equation models. They 
will be described here briefly. The algebraic models are the 
Thomas, Baldwin-Lomax, and Modified Mixing Length- 
Thomas (MMLT) models. The Thomas model (based on 
the work of Ref. 7) is the standard algebraic turbulence 
model in PARC. This model calculates turbulent viscosity 
near surfaces (wall-bounded part of the model) and in 
regions where flows are mixing (free-shear layer part of the 
model) but is optimized for the latter. The Baldwin-Lomax 
algebraic turbulence model 8 is also available in PARC. This 
model only calculates turbulent viscosity in wall-bounded 
regions. The third algebraic model uses the Modified 
Mixing Length (MML) model (originally developed to 
analyze iced airfoils) 9 for wall-bounded regions and the 
Thomas model for free-shear layer regions. This combina- 
tion turbulence model was developed from its two compo- 
nents in Ref. 10 and is referred to as the MMLT model for 
the rest of this discussion. 

Algebraic turbulence models such as those described 
here often model complex flows inadequately because they 
use single mixing length distributions to calculate turbulent 
viscosity, which often are not applicable to all flows. Two- 
equation models avoid this single mixing length limitation 
by using additional transport equations to calculate turbulent 
viscosity. However, these models are substantially more 
computationally expensive. The two-equation models that 
have been installed in PARC and were investigated in this 
study are the Chien low Reynolds number k-£ model 11 with 
modifications for compressibility added by Nichols 12 and the 
low Reynolds number k-£ model based upon the work of 
Speziale. 13 

Discussion of Flow Cases 

The transonic diffuser flows considered in this study 
are those described in Refs. 14-17. This two-dimensional 
diffuser with a convergent-divergent channel was designed 


to simulate the types of flows that exist in supersonic inlets 
of aircraft engines. Extensive flow measurements were 
made during tests of this diffuser for flows with and without 
externally applied oscillations. Only the experimental 
results for unexcited flows were examined for comparison 
with PARC flow calculations. 

A schematic of the diffuser geometry is shown in 
Fig. 1 . The diffuser had an entrance-to- throat area ratio of 
1 .4 and an exit-to- throat area ratio of 1 .5. The distance 
between sidewalls was approximately four throat heights. 
Suction slots were placed in the sidewalls and top comers to 
keep the flow two-dimensional. Three flows were investi- 
gated with PARC2D using the five turbulence models. 

These flows were defined by the ratio of the exit static 
pressure to the inflow total pressure (R) and are referred to 
as the no-shock (R = 0.862), weak-shock (R = 0.82), and 
strong-shock (R = 0.72) cases. The inflow total temperature 
was approximately 300 K and the outflow static pressure 
was atmospheric for all cases. 

No-Shock Case 

The first flow examined with PARC was that with no 
shock forming in the duct. The back pressure (R = 0.862) 
was high enough to prevent supersonic flow from forming 
in the diffuser downstream of the throat. A grid sensitivity 
investigation was conducted using this no-shock case before 
PARC solutions using the five turbulence models were 
compared with each other and with experimental data. This 
grid sensitivity study is discussed next. 

Similar numerical studies 18 * 22 used a computational 
grid having 81 points in the horizontal direction and 51 
points in the vertical direction (81 x 51). The two grids 
constructed for this investigation are shown in Fig. 2. The 
first grid had 81x51 points with the point next to either 
wall placed in the laminar sublayer (y+ < 5). The second 
grid also had 81 points in the horizontal direction but a total 
of 81 points in the vertical direction. These extra points in 
the second grid were used to pack the boundary layer 
regions more tightly with the point next to the wall at a 
distance corresponding to y ~ 1. For each turbulence 
model, the solution obtained with the 81 x 51 grid was 
compared to the solution obtained with the 81 x 81 grid to 
determine grid sensitivity. 

Only the Speziale k-£ model results showed signifi- 
cant differences between solutions obtained with the two 
grids. Low Reynolds number k-£ turbulence models often 
require more tightly packed grids than other turbulence 
models, such as the three algebraic models investigated in 
this study. Not only should the first grid point from the wall 
be placed in the laminar sublayer, it should correspond to a 
position of y+ ~ 1. Avva, Smith, and Singhal 23 report that if 
the first grid point is placed in the logarithmic layer instead 
of the laminar sublayer when using low Reynolds number 
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k-E turbulence models, the peak value of turbulent kinetic 
energy can increase by a factor of 2, thus creating inaccurate 
flow predictions. Figure 3 shows a comparison of turbulent 
kinetic energy profiles obtained from the Chien and Speziale 
k-£ solutions along the bottom wall of the diffuser at the 
throat. The Chien model produces similar profiles (see Fig. 
3a) for solutions obtained with the two grids. The turbulent 
viscosity was much lower than the molecular viscosity at the 
first point off the wall for both grids. This indicates that the 
first points were in the laminar sublayer. 

The Speziale model results (Fig. 3b) show that the 
turbulent kinetic energy in the near wall region calculated 
with the 81 x 51 grid is substantially higher than that 
calculated using the 81 x 81 grid. Also, the turbulent 
viscosity at the first point off the wall was higher than the 
molecular viscosity for the 81 x 51 grid, indicating that the 
first point was outside the laminar sublayer. The turbulent 
viscosity calculated with the 81 x 81 grid did not exceed the 
molecular viscosity until the third or fourth point away from 
the wall. Based upon this grid sensitivity investigation, the 
81x51 grid is sufficient for the three algebraic turbulence 
models and the Chien k-E model, but the 81 x 81 grid is 
required for the Speziale k-E model to obtain calculations 
that compare with the experimental data. A comparison of 
PARC solutions with experimental data follows next. 

Pressure distributions obtained from the PARC 
calculations are compared with experimental pressures 
along the top and bottom walls of the diffuser in Fig. 4. The 
Baldwin-Lomax and Thomas solutions predicted pressures 
near the throat (x/Rj. ~ 0) that are too low compared with the 
experimental data whereas the Speziale k-E solutions 
predicted pressures that are too high. The Chien and 
MMLT solutions predicted pressures in this location more 
accurately than the other models. Mach number gray-scale 
contours for the no-shock flow (Chien k-E solution) are 
shown in Fig. 5 and the maximum Mach numbers in the 
diffuser using each turbulence model are listed in Table 1 . 
These maximum Mach numbers (occurring at the throat) 
correspond to the pressure distributions of Fig. 4. That is, 
the Baldwin-Lomax and Thomas solutions have the highest 
maximum Mach numbers corresponding to their lowest 
pressures at the throat. Also, the Speziale solution produces 
the lowest maximum Mach number in the diffuser corre- 
sponding to its highest pressures at the throat. 


Table 1 Calculated flow properties for the no-shock flow. 


Case 

Maximum Mach 
number in duct 

Thomas 

0.863 

Baldwin-Lomax 

.872 

MMLT 

.828 

Chien k-e 

.813 

Speziale k-e 

.779 


Convergence was determined when the residual error 
dropped several orders of magnitude and when the solutions 
did not change with more iterations. Convergence histories 
for all the solutions are listed in Table 2. Each solution took 
10 000 iterations to converge, but the k-E solutions required 
more CPU time than the algebraic model solutions. The two 
k-E model solutions took the same amount of CPU time to 
obtain the same number of iterations even though the 
Speziale calculations required a larger grid (81x81) than 
the Chien calculations (81x51). For the implementations 
used in PARC, the Speziale model requires less CPU time 
per iteration per grid point than the Chien model. 


Table 2.- Convergence histories for the no-shock flow. 


Case 

Iterations 

Cray Y/MP 
CPU time (s) 

Thomas 

10 000 

600 

Baldwin-Lomax 

10 000 

600 

MMLT 

10 000 

600 

Chien k-e 

10 000 

1700 

Speziale k-e 

10 000 

1700 


Weak-Shock Case 

The next flow examined was that producing a weak 
shock in the diffuser. This flow was initialized by setting 
the outflow static pressure to a very small value ( R = 0.12) 
to allow supersonic flow to exist from the throat to the exit 
plane. The back pressure was then set to its correct level for 
this flow (R = 0.82) to allow a shock to form slightly 
downstream of the throat (see Mach number contours for the 
Chien solution in Fig. 6). As with the no-shock case, these 
weak-shock calculations were continued until the solutions 
did not change. The experimental results described in Refs. 
14 and 15 indicated some very small self-induced oscilla- 
tions about the mean flow. All the PARC calculations 
reached non-oscillatory solutions, as expected since PARC 
was run in steady-state mode. 

Pressure distributions along the top and bottom walls 
of the diffuser from the PARC calculations are compared 
with the mean flow experimental pressures in Fig. 7. All 
calculated pressure distributions agree relatively well with 
the experimental data except for the Speziale solution, 
which predicts the pressure to rise on the top and bottom 
walls much earlier than the other solutions or experimental 
data. 


Velocity profiles obtained from the PARC calcula- 
tions are compared to laser Doppler velocimeter data 
(documented in Ref. 14 for the weak-shock case discussed 
here and the strong-shock case discussed in the next section) 
in Fig. 8. At x/Rj. = 1.729, the Baldwin-Lomax solution 
indicates that the flow is just passing through the weak 
shock while the other solutions and data indicate that the 
shock occurred upstream. As a result, the Baldwin-Lomax 
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velocities are nearly sonic in the core of the flow at this 
location (Fig. 8a) while the core flow velocities of the other 
solutions are subsonic. 

The maximum Mach numbers in the diffuser for each 
solution are listed in Table 3. The experimental value listed 
is the Mach number just outside the upper- wall boundary 
layer. The positions of the shock in the core of the flow 
obtained from the PARC calculations and the experimental 
data are also given in Table 3. For the PARC calculations, 
the positions of the shocks were determined to be the 
locations where the maximum drop in Mach number from 
any streamwise point to the next downstream point oc- 
curred. These shock positions correspond to the pressure 
distributions of Fig. 7 and the velocity profiles of Fig. 8. 

The Speziale solution, for example, predicted both the shock 
and the associated rise in static pressure to occur upstream 
of the other PARC solutions and experimental data while 
the Baldwin-Lomax solution predicted both to occur further 
downstream of the other solutions and experimental data. 


Table 3.- Calculated flow properties for the weak-shock 
flow. 


Case 

Shock position 
(throat heights, 

Maximum Mach 
number in duct 

Thomas 

1.615 

1.285 

Baldwin-Lomax 

1.690 

1.298 

MMLT 

1.537 

1.255 

Chien k-£ 

1.283 

1.216 

Speziale k-e 

.884 

1.153 

Salmon Data 

1.47 

1.235 (upper wall) 


Convergence information from the PARC calculations 
is presented in Table 4. The Baldwin-Lomax solution took 
the most iterations to achieve steady state. The k-e solutions 
required the fewest iterations, but the most CPU time. 


Table 4.- Convergence histories for the weak-shock flow. 


Case 

Iterations 

Cray Y/MP 
CPU time (s) 

Thomas 

25 000 

1400 

Baldwin-Lomax 

30 000 

1800 

MMLT 

25 000 

1500 

Chien k-e 

16 000 

2800 

Speziale k-e 

16 000 

2800 


Strong-Shock Case 


The last diffuser flow examined with PARC was the 
flow forming a strong shock (R = 0.72). These calculations 
were also begun with a very low back pressure (R = 0.12) to 
set supersonic flow in the diffuser from the throat to the exit 
plane. The back pressure was then increased to the correct 
outflow static pressure. After each calculation was run for 
30 000 iterations, only the Chien and Speziale calculations 
reached steady solutions. The three calculations using 
algebraic turbulence models demonstrated oscillations in the 
flow that did not decay with more iterations, although 
PARC was still being run in steady-state mode. 


When it was determined that the algebraic turbulence 
model solutions were not going to reach steady state, the 
calculations were started again using the same initial 
conditions with supersonic flow from the throat to the exit 
plane. All 3 calculations obtained with the algebraic 
turbulence models were run for 8000 iterations, which was 
sufficient to allow the strong shock to form in the general 
region of the diffuser where the shock was oscillating during 
the previous iterations. At this point, all three calculations 
were run for 8 sets of 4000 iterations each (32 000 iterations 
for each calculation after the initial 8000 iterations) to obtain 
8 intermediate solutions for each calculation that could be 
averaged to provide mean flow properties. References 14 
and 1 8 indicate that the magnitudes of the self-induced 
oscillations for this strong-shock case are considerably 
higher than those for the weak-shock case. This natural 
unsteadiness in the flow is probably a major reason for the 
inability of PARC (using the algebraic turbulence models) 
to produce steady solutions. 


Figure 9 compares the PARC results to the experi- 
mental mean flow pressure distributions. Pressure distribu- 
tions for each algebraic turbulence model were obtained by 
averaging the eight intermediate solutions stored for each 
model. The k-e solutions (the only solutions not averaged 
because they were not oscillatory) match the experimental 
pressure profiles best. The experimental data indicate that 
the flow separates from the top wall because of the strong 
shock and reattaches at about six throat heights downstream 
of the throat. All the PARC solutions predict the flow 
separation from the top wall due to the shock, but only the 
k-e solutions predict reattachment within the computational 
domain. This may be observed in the velocity profiles 
shown in Fig. 10. The Chien solution predicts a reattach- 
ment location at 5.4 throat heights downstream of the throat 
and the Speziale solution predicts a reattachment at 5.9 
throat heights. At x/H^ =6.340 (Fig. 10c)andx/H r = 
7.493 (Fig. lOd) the algebraic solutions all have large flow 
reversal regions near the top wall. These much larger flow 
separations also contributed to the inability of the algebraic 
solutions to converge. 


Figure 1 1 shows gray-scale Mach number contours 
for the Chien solution and Table 5 gives the maximum 
Mach numbers in the duct and the shock positions for the 
PARC solutions and the experimental data. For the 
algebraic turbulence model solutions, the shock positions 
are averaged. The standard deviations of the pressure 
distributions for the algebraic model solutions are shown in 
Fig. 12. Overall, the Baldwin-Lomax solution shows the 
largest variation, particularly near the shock location, while 
the Thomas solution has the smallest variation. The 
convergence histories for the solutions are given in Table 6. 
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Table 5.- Calculated flow properties for the strong-shock 
flow. 


Case 

Shock position 
(throat heights, 

Maximum Mach 
number in duct 

Thomas 

2.413 

1.422 

Baldwin-Lomax 

2.224 

1.383 

MMLT 

2.350 

1.405 

Chien k-e 

2.444 

1.416 

Speziale k-e 

2.161 

1.373 

Salmon Data 

2.39 

1.353 (upper wall) 


Table 6.- Convergence histories for the strong-shock flow. 


Case 

Iterations 

Cray Y/MP 
CPU time (s) 

Thomas 

40 000 

2400 

Baldwin-Lomax 

40 000 

2400 

MMLT 

40 000 

2400 

Chien k-e 

30 000 

5000 

Speziale k-e 

30 000 

5000 


Conclusions 

An evaluation of five turbulence models available in 
the PARC code, an internal flow Navier-Stokes solver used 
extensively by the propulsion community, has been con- 
ducted for flows through a two-dimensional transonic 
diffuser. These models are the Thomas, Baldwin-Lomax, 
and Modified Mixing Length Thomas (MMLT) algebraic 
turbulence models and the Chien and Speziale low Reynolds 
number k-e (two-equation) models. Three diffuser flows, 
having different flow conditions characterized by the ratio 
of the outflow static pressure to the inflow total pressure, 
were calculated. These flows, which are representative of 
many inlet cases to which PARC has been applied, were 
referred to as the no-shock, weak-shock, and strong-shock 
cases. 

Overall, the Chien k-e model was the most accurate of 
the five turbulence models when all three flows are consid- 
ered. However, the MMLT model provided results as 
accurate as those of the Chien model for the no-shock case 
and the weak-shock case at a significantly lower computa- 
tional cost (measured in CPU seconds required to provide 
converged solutions). For the strong-shock case, only the 
Chien and Speziale k-s models produced steady-state 
solutions. None of the PARC calculations using the 
algebraic turbulence models reached steady state, but 8 
intermediate solutions were obtained at 4000 iteration 
increments for each model and averaged to provide mean 
flow properties. The experimental data indicates that the 
strong-shock flow demonstrates large scale unsteadiness, 
unlike the no-shock and weak-shock flows. This probably 
contributes to the inability of the algebraic model calcula- 
tions to reach steady state. The comparison of the averaged 
algebraic model solutions and steady k-e solutions to the 
experimental mean flow properties indicates that the two k-e 
solutions matched the experimental data best. 


As improvements in turbulence models continue to be 
made, the availability of several turbulence models (e.g. 
both algebraic and two-equation models) in a multipurpose 
Navier-Stokes code like PARC allows selection of the 
turbulence model appropriate for the flow to be analyzed. 
For flows where an algebraic model may be expected to 
provide results as accurate as a more sophisticated turbu- 
lence model (like k-e), using the algebraic turbulence model 
would save computational resources and require less time to 
obtain solutions. For more difficult flows where an 
algebraic model would not produce flow predictions as 
accurate as a more complex model, using the more complex 
model would be justified. 
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Figure 1 - Diffuser geometry. 



b. 81 x 81 grid. 

Figure 2 - Computational grids for diffuser flow cases. 



Figure 3 - Turbulent kinetic energy profiles at the throat (no-shock case). 
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Figure 4 - Pressure distributions for the no-shock case. 


Mach Number 
Contours 



Figure 5 - Mach number contours for the no-shock case (Chien k-e solution). 


Mach Number 
Contours 



Figure 6 - Mach number contours for the weak-shock case (Chien k-£ solution). 
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Figure 7 - Pressure distributions for the weak-shock case. 




a. x/H = 1.729 


b. x/H = 2.882 




c. 


Figure 8 - Velocity profiles for the weak-shock case. 
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Figure 9- Pressure distributions for the strong-shock case. 






c. x/H = 6340 d. x/H = 7.493 

Figure 10 - Velocity profiles for the strong-shock case. 
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d/d 



Figure 11 - Mach number contours for the strong-shock case (Chien k-e solution). 




Figure 12- Standard deviations of pressure for the strong-shock case. 
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